Cosmic-Ray Rejection by Laplacian Edge Detection 
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ABSTRACT 

Conventional algorithms for rejecting cosmic-rays in single CCD exposures rely on the con- 
trast between cosmic-rays and their surroundings, and may produce erroneous results if the Point 
Spread Function (PSF) is smaller than the largest cosmic-rays. This paper describes a robust al- 
gorithm for cosmic-ray rejection, based on a variation of Laplacian edge detection. The algorithm 
identifies cosmic-rays of arbitrary shapes and sizes by the sharpness of their edges, and reliably 
discriminates between poorly sampled point sources and cosmic-rays. Examples of its perfor- 
mance are given for spectroscopic and imaging data, including Hubble Space Telescope WFPC2 
images. 

Subject headings: instrumentation: detectors — methods: data analysis — techniques: image processing 



1. Introduction 

Various methods are in use for identifying and 
replacing cosmic-ray hits in CCD data. The 
most straightforward approach is to obtain mul- 
tiple exposures of the same field (or multiple non- 
destructive readouts during a single exposure; e.g., 
Fixsen et al. 2000). In general, a given pixel 
will suffer from a cosmic-ray hit in only one or 
two of the exposures, and remaining exposures 
can be used to obtain its replacement value (e.g., 
Zhang 1995). Methods for combining multiple ex- 
posures have reached a high degree of sophisti- 
cation, particularly those developed for dithered 
Hubble Space Telescope (HST) data (e.g., Wind- 
horst, Franklin, & Neuschaefer 1994, Freudling 
1995, Fruchter & Hook 1997). 

However, there are circumstances when cosmic- 
ray identification in single exposures is required 
or desirable. The object of interest may be vary- 
ing or moving on short timescales, and in the case 
of long-slit spectra the positions and intensities 
of sky lines and object spectra may change (e.g., 
Croke 1995). Furthermore, pixels can be hit by 
cosmic-rays in more than one exposure, and some 
affected pixels may remain after combining indi- 
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vidual images. Cosmic-ray removal in individual 
exposures may also be desirable if the images are 
shifted with respect to each other by a non-integer 
number of pixels, or if the seeing varied signifi- 
cantly between the exposures (see Rhoads 2000). 
Finally, multiple exposures are simply not always 
available. 

Methods for identifying cosmic-rays in single 
images or spectra include median filtering (e.g., 
QzAP by M. Dickinson), filtering by adapted Point 
Spread Functions (PSFs) (e.g., Rhoads 2000), 
trained neural networks (Salzberg et al. 1995), 
and interpolation of neighbouring pixels (e.g., the 
COSMICRAYS task in the IRAF package). AU these 
methods effectively remove small cosmic-rays from 
well sampled data. 

The most widely used methods are based on 
some form of median filtering, and usually include 
adaptations to exclude stars and emission lines 
from the list of cosmic-rays. However, problems 
arise when cosmic-rays affect more than half the 
area of the filter, or the PSF is smaller than the 
filter. The size of the filter is therefore a trade- 
off between detecting large cosmic-rays and limit- 
ing contamination by structure on the scale of the 
PSF. 

In this paper, a new algorithm for rejecting 
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cosmic-rays in single exposures is described. It 
is based on Laplacian edge detection, which is a 
widely used method for highlighting boundaries in 
digital images (see, e.g., Gonzalez & Woods 1992). 
The strength of the method is that it relies on the 
sharpness of the edges of cosmic-rays rather than 
the contrast between entire cosmic-rays and their 
surroundings. Therefore, it is largely independent 
of the morphology of cosmic-rays. This property 
is very useful, and forms the basis for a robust dis- 
crimination between poorly sampled point sources 
and cosmic-rays. 

2. The Laplacian 

The Laplacian of a 2-D function is a second- 
order derivative defined as 

The Laplacian is commonly used for edge detec- 
tion in digital images. In this application the im- 
age is convolved with the Laplacian of a 2-D Gaus- 
sian function of the form 

f{x,y) = e^v{-^y (2) 

where = x^+y"^ and ct is the standard deviation. 
The second-order derivative with respect to r has 
the form 
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The Laplacian has zero-crossings at r = ±-\/2cr, 
and the locations of edges are found by identifying 
zero-crossings in the convolved image. The stan- 
dard deviation can be tuned to the smoothness of 
the edge, and by using a range of values for a both 
sharp and smooth edges can be identified (Marr & 
Hildreth 1980). 

Cosmic-rays have very sharp edges, and the 
convolution kernel should be most sensitive to 
variations on small scales. The appropriate dis- 
crete implementation of Eq. 3 has the form 



(4) 



The average value of a Laplacian image (obtained 
by convolving an image with the kernel given in 
Eq. 4) is zero, and smooth structure in the image 
is removed. 



3. Implementation 
3.1. Basic Procedure 

A straightforward convolution of Eq. 4 with 

a CCD image produces negative cross patterns 
around high pixels. As a result, connected cosmic- 
ray pixels suffer from attenuation by the negative 
cross patterns of their neighbors. Hence before 
convolution the original image / needs to be sub- 
sampled. The results are independent of the sub- 
sampling factor, and a factor two is computation- 
ally least expensive: 
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with nxn the size of the original image and i, j = 
1, . . . , 2n. The Laplacian of the subsampled image 
is 

£(2) = vVo-r^'\ (6) 

with o denoting convolution. 

The Laplacian of the edge of a cosmic-ray is 
negative on the outside and positive on the inside 
of the cosmic-ray. Therefore, by setting all nega- 
tive values in the Laplacian image to zero cosmic- 
ray affected pixels are retained and their negative 
cross patterns are removed: 
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Finally, the image is resampled to its original res- 
olution: 
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with i, j — I, . . . , n. 

The process is illustrated in Fig. 1. A small 
section of a 2D long slit spectrum is shown in (a). 
Panel (b) shows the same image, subsampled by a 
factor six and convolved with the Laplacian kernel 
given in Eq. 3. Edges are at the locations of zero- 
crossings. In panel (c) all negative values are set 
to zero. Finally, the image is block averaged by 
a factor six (panel d). The cosmic-ray stands out 
clearly. Because the spectrum is smooth on scales 
of ~ 1 pixel its Laplacian is close to zero. 

The numerical value associated with each edge 
is the difference between the two neighbouring pix- 
els, and the resampled Laplacian of an image / 
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Fig. 1. — Illustration of Laplacian edge detection. 
The original image is shown in (a). Panel (b) shows 
the same image after subsampling by a factor six and 
convolution with the Laplacian kernel. Edges are pos- 
itive on the inside of the cosmic-ray, and negative on 
the outside. Negative pixels are set to zero in (c), and 
the image is block averaged to its original resolution 
in (d). 
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Fig. 2. — The thick solid line shows the noise proper- 
ties of a simulated image with ~ 10^ 7tT cosmic-rays. 

The hatched histogram shows the distribution of pixel 
values in the Laplacian image. Pixel values in the 
Laplacian image were divided by the subsampling fac- 
tor. In the positive tail of the noise distribution the 
characteristics of the original image are (to a good ap- 
proximation) conserved. The width of the distribution 
of cosmic-rays is increased by 0.5(7. 



consisting of a smooth background B with super- 
posed noise and cosmic-rays is approximately 
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with fg the subsampling factor that was used. The 
Laplacian thus retains the flux in high pixels and 
removes the local background, making it very use- 
ful for identifying cosmic-rays. 

To identiiy cosmic-rays in the Laplacian image 
the value of each pixel is compared to the expected 
noise at that location. The noise characteristics of 
the Laplacian image arc not the same as in the 
original image: the Laplacian operator itself in- 
creases the noise by a factor y^5/4, and all nega- 
tive Poisson fluctuations in the original image arc 
(close to) zero in the Laplacian image. These ef- 
fects are demonstrated in Fig. 2 for a simulated 
image containing ^ 10^ 7a cosmic-rays. As ex- 
pected, the distribution of pixel values near zero 
is strongly distorted by the Laplacian. However, 
the positive tail of the distribution (i.e., the region 



where cosmic-rays are found) remains virtually un- 
changed. Therefore the noise properties of the 
original image can be used for identifying cosmic- 
rays in £+, which greatly simplifies the analysis. 

The noise model is constructed by convolving 
the original image with a median filter: 



A^ = 5"V5(M5o7)+<72„ 



(10) 



where g is the gain in electrons/ ADU, am is the 
readnoise in electrons, and M5 is a 5 x 5 median 
filter. For long slit spectroscopic data the algo- 
rithm offers the option of fitting and subtracting 
sky lines and/or object spectra. These fits provide 
the basis for the long slit noise model, which is op- 
timized by applying Eq. 10 to the residuals from 
the fits. 

The Laplacian image is divided by the noise 
model to obtain the deviations from the expected 
Poisson fluctuations: 
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Cosmic-rays arc identified by selecting pixels in S 
that are above a given threshold ciim- 

For single-pixel cosmic-rays on a smooth back- 
ground the detection probability depends on the 
noise in the background (Torg only. The error in the 
background estimate scales clS (Tqyp: I yjn, with n the 
number of pixels. Since four neighbouring pixels 
are used to determine the Laplacian of a given 
pixel, the width of the distribution of cosmic-rays 
is increased by crorg/2 (see Fig. 2). As an example, 
if a detection threshold of 5cr is applied, ~ 5 % of 
4cr peaks in the original image will be marked as 
cosmic-rays, ~ 50% of 5cr peaks, and 95% of 
6(7 peaks. 

The detection probability of cosmic-rays larger 
than a single pixel depends on the number of 
connected pixels and the pixel-to-pixel variation 
within the cosmic-ray. In general the Laplacian 
(and hence the detection probability) is lowest for 
cosmic-rays with small pixel-to-pixel variations. 
In the limiting case of a cosmic-ray with negligible 
pixel-to-pixel variation 

S^,i-N-](}-'^){h,-B,^,\ (12) 

with n^j- the number of cosmic-ray pixels adja- 
cent to pixel (i, j). Pixels on the corners of large 
cosmic-rays have at most two adjacent pixels, and 
the Laplacian always retains at least « 50 % of 
their flux. Therefore, (arbitrarily) large cosmic- 
rays can be removed by applying the rejection 
process iterativcly. Before each iteration, pixels 
flagged as cosmic-rays in previous iterations are 
replaced by the median of surrounding "good" pix- 
els, acd 

3.2. Removal of Sampling Flux 

The Laplacian gives the difference between one 
pixel and its neighbours, and contains no informa- 
tion on the nature of detected features. Real as- 
tronomical objects produce signal in the Laplacian 
image because of Poisson noise and because their 
intrinsically smooth intensity profiles are sampled 
by the pixels. This "sampling flux" is generally 
small, and will not produce spurious detections as 
long as it does not exceed the predicted Poisson 
fluctuations (see Eq. 10, 11). However, for bright 
objects the sampling flux can be significant, in par- 
ticular if the Point Spread Function (PSF) is not 
well sampled. 



Sampling flux is removed from S in two steps. 
First, all structure that is smooth on scales of > 5 
pixels is removed by a 5 x 5 median filter: 

5' = ^- (50M5). (13) 

This procedure very effectively removes sampling 
flux resulting from extended bright objects (in- 
cluding point sources if the PSF is well sampled). 
Because of the large size of the filter cosmic-rays 
and the noise properties of S remain unaffected. 

Next, sampling flux resulting from critically 
sampled (or even undersampled) point sources is 
removed. As is well known, it is very hard to 
distinguish cosmic-rays from stars and emission 
lines in marginally sampled data, because they can 
have very similar pixel-to-pixel variations within 
an area < 3 x 3 pixels. 

Point sources are distinguished from cosmic- 
rays by their symmetry. An image only containing 
symmetric fine structure on scales of 2 — 3 pixels is 
constructed and compared to the Laplacian image. 
This "fine structure image" T is created from the 
original image by a combination of median filters: 

jr=(M3o7)-([M3o7]oM7), (14) 

where M„ is an n x n median filter. The second 
term serves to remove large scale structure from 
T. An important property of T is that central 
pixels of undersampled point sources do not van- 
ish, but are replaced by the median of the sur- 
rounding pixels. The Laplacian image is divided 
by the fine structure image, and cosmic-rays are 
selected as those pixels which have S' > aum and 
L'^ jT > /lim, with /lim defining the minimum 
contrast between the Laplacian image and the fine 
structure image. 

Figure 3 demonstrates the procedure. Artificial 
images of three stars and a cosmic-ray arc shown 
in the top row. One stellar image is slightly over- 
sampled (cr = 1.5 pixels), the second is critically 
sampled (cr = 1.0 pixels), and the third is slightly 
undersampled (tr = 0.7 pixels). The contrast be- 
tween the cosmic-ray and its local background is 
identical to the contrast between the central pixel 
of the undersampled star and its local background. 
As a result, the highest pixels in their Laplacian 
images (shown in the second row) have very 
similar values. The third row shows fine struc- 
ture images J^. The cosmic-ray has very low sig- 
nal in J^, whereas the stars retain a significant 
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fraction of their flux because of their symmetry. 
The bottom row shows the Laplacian images di- 
vided by the fine structure, C^' jT. Only pixels 
with / T > /lini arc retained. The critically 
sampled star has L'^ jT = 0.7, and L'^ jT = 1.8 
for the undersampled star. The cosmic-ray has 
£+/jF = 21, and is easily distinguished from un- 
dersampled point sources. 
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Fig. 3. — Differentiating between marginally sampled 
point sources and cosmic-rays. Panels show, from top 

to bottom, artifical images of stars and a cosmic-ray, 
the Laplacian of these images , their fine structure 
image T, and the Laplacian divided by the fine struc- 
ture /T. The number in each panel is the value of 
the highest pixel. The highest pixels in the Laplacian 
images of the undersampled star (a = 0.7 pixels) and 
the cosmic-ray are similar. However, they are very 
different after division by the fine structure image. 

In general, the appropriate value of /lim de- 
pends on the sampling of a given point source, 
its S/N ratio, and whether it lies on the center 
of a pixel or close to the edge. These effects 
can be simulated by creating artificial PSFs with 
varying sampling, S/N ratios, and subpixcl posi- 
tions, and calculating £+/jF for each star. Fig- 
ure 4 shows the dependence of jT on the sam- 
pling, expressed as the hdl width at half maxi- 
mum (FWHM) of stars (in pixels). The width of 
the shaded region demonstrates the (maximum) 
effects of varying S/N ratio and subpixel position. 



These simulations demonstrate that the value of 
L'^ jT , and hence the appropriate choice of /um, 
depends mainly on the sampling. For data that 
arc well sampled /lim = 2 is appropriate; hence 
this is the default value in the algorithm. For un- 
dersampled data higher values of /nm are needed 
to discriminate point sources and cosmic-rays; as 
an example, data taken with the Wide Field chips 
in HST's WFPC2 camera require /lim ~ 5. 
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Fig. 4. — Dependence of the ratio of the Laplacian 
and the fine structure image T on the sampling, 
as derived from simulations (see text). For a given 
FWHM of stars the value of /um should be chosen 
such that it exceeds the dashed band. The default 
/lim = 2 is appropriate for data that are critically, or 
better, sampled. HST WFPC2 data require /um « 5. 



3.3. Additional Features 

The basic algorithm detects cosmic-rays and re- 
jects point sources from the initial list. In ad- 
dition, the program L.A.COSMIC allows a lower 
detection threshold to be used for pixels neigh- 
bouring those already flagged as cosmic-rays. It 
also replaces cosmic-rays by the median of sur- 
rounding "good" pixels, and offers the option of 
applying the algorithm iteratively. On a Sun Ul- 
traSparc 1 (200 MHz) the IRAF implementation of 
L. A. Cosmic requires approximately 65 s per iter- 
ation for an image of 800 x 800 pixels. The run 
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time scales linearly with the number of pixels. 

4. Examples 

The algorithm was tested on a variety of real 
and artificial data sets, consistently producing 
very good results. The examples given here serve 
to illustrate its performance. 

4.1. Well Sampled Imaging CCD Data 

Figure 5(a) shows a well sampled artifical im- 
age containing 500 stars, 100 galaxies, and 227 
cosmic-rays. Stars have a — 1.5 pixels, equivalent 
to, for example, FWHM = 0'.'78 seeing with 0'.'22 
pixels. All cosmic-rays arc > 5(T above the sky 
background. The reconstruction by L.A.COSMIC 
is shown in (b). Panel (c) shows the input cosmic- 
ray image, and panel (d) shows the cosmic-rays 
found by L.A.COSMIC. The program found 222 
of the 227 cosmic-rays (98%). Importantly, only 
1 of the 500 stars (0.2 %) and none of the galaxies 
was inadvertently identified as a cosmic-ray. 

For well sampled imaging data the performance 
is similar to median filtering methods such as 
QzAP (by M. Dickinson). Sophisticated median 
filtering recovers close to 100 % of cosmic-rays that 
are smaller than the filter size, but breaks down if 
the cosmic-ray is larger than the filter, and/or the 
FWHM of the PSF is smaller than the filter. 

4.2. HST WFPC2 Data 

Cosmic-rays in images obtained with WFPC2 
on HST are notoriously difiicult to remove, be- 
cause of their large number and the undersampling 
of the PSF. Nevertheless, the algorithm described 
here performs very well on WFPC2 data. The 
method is insensitive to the size of cosmic-rays, 
and the undersampling of the PSF can be taken 
into account by setting the parameter /um = 5 
(see Sect. 3.2). 

Figure 6(a) shows part of a 2400 s WF obser- 
vation in /f814W of galaxy cluster MS 1137+67. 
The reconstruction of the image by L.A.COSMIC 
is shown in (b). 



Virtually all cosmic-rays are removed, and none 
of the real objects is mistaken for a cosmic-ray. 
The small panels show examples of stars and 
galaxies in WF chips, extracted from WFPC2 ob- 
servations of various targets. The algorithm leaves 
stars intact, and is able to remove arbitrarily large 
cosmic-rays. 

The reliability of cosmic-ray identification can 
be tested by comparing the results of L.A.COSMIC 
on single images to "true" cosmic-ray images cre- 
ated from multiple exposures. The test used one 
of the CR-SPLIT WFPC2 images of the cluster 
MS 2053-04 {z = 0.58). As a result of its low 
Galactic latitude approximately 50 % of /f814w < 
22 objects are stars. The reduction of these data 
is described in van Dokkum et al. (2001). The 
algorithm found 5638 (98.1 %) of 5750 cosmic-ray 
affected pixels deviating more than 6cr from the 
background, and 4687 (99.1%) of 4729 pixels de- 
viating more than IOct. The number of false pos- 
itives, i.e., pixels inadvertently marked as cosmic- 
rays, is 72 (1.2%) at > 6(7, and only 1 (0.02%) 
at > 10(7. These numbers compare favorably to 
cosmic-ray rejection algorithms based on morpho- 
logical classification by neural networks (Salzberg 
et al. 1995). 

4.3. Spectroscopic CCD Data 

The algorithm optimized for spectroscopic long- 
slit data is very similar to the implementation for 
imaging data. The main difference is that the 
program offers the possibility of fitting and sub- 
tracting sky lines and the object spectrum before 
convolution with the Laplacian kernel. 

An example (an 1800 s long slit spcctnun of a 
galaxy, obtained with the Low Resolution Imag- 
ing Spectrograph on the W. M. Keck Telescope) 
is shown in Fig. 7. Note that the Laplacian is 
very effective in removing the fringe pattern that 
is present after subtracting strong sky lines. The 
fine structure image T is used to identify emission 
lines and other sharp features in the spectra, in 
similar fashion as the identification of undersam- 
pled stars in imaging data. 
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Fig. 5. — (a) Artificial image containing 500 stars, 100 galaxies, and 227 cosmic-rays, (b) Reconstruction of the 
image by L. A. Cosmic. The true cosmic-ray image is shown in (c), and the cosmic-rays found by L. A. Cosmic are 
shown in (d). 
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Fig. 6. — (a) HST WFPC2 image of galaxy cluster MS 1137+67. The restoration by L. A. Cosmic is shown in (b). 
Small panels show close-ups for a selection of stars and galaxies in various WFPC2 images. The algorithm leaves 
stars intact, and is able to remove cosmic-rays of arbitrary shapes and sizes. 
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Fig. 7. — Demonstration of Laplacian cosmic-ray rejection for long-slit spectra. Prom top to bottom are shown: the 
original spectrum, the residuals after subtraction of ID-fits to the sky lines and the object spectrum, the Laplacian 
image, the Laplacian image divided by the noise model, the finestructure image, the bad pixel map, and the restored 
spectrum. Note that the Laplacian very effectively removes fringing. In this case all cosmic-rays are removed in a 
single iteration. The bright emission line is not marked as a cosmic-ray because of its prominence in the fine structure 
image. 
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5. Conclusions 

Cosmic-rays in single images or spectra can be 
removed by a variation of Laplacian edge detec- 
tion. The procedure is robust, and requires very 
few user-defined parameters. The method rejects 
cosmic-rays of arbitrary size and distinguishes un- 
dersampled point sources from cosmic-rays with 
high confidence. It is implemented in the pro- 
gram L. A. Cosmic, which can be obtained from 
http://www.astro.caltech.edu/~pgd/lacosmic/ . 
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